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Abstract 

This paper is concerned with the reduction of a noisy synchronous 
Boolean network to a coarse-grained Markov chain model. Consider an n- 
node Boolean network having at least two basins of attraction and where 
each node may be perturbed with probability < p < 1 during one time 
step. This process is a discrete-time homogeneous Markov chain with 2" 
possible states. Now under certain conditions, the transitions between the 
basins of the network may be approximated by a homogeneous Markov 
chain where each state of the chain represents a basin of the network, i.e. 
the size of the reduced chain is the number of attractors of the original 
network. 
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1 Noisy Boolean networks (NBNs) 
1.1 Definition 

Boolean networks (BNs) have been used for several decades as models of bio- 
chemical networks, mainly to predict their qualitative properties and, in the 
case of genetic regulatory network s, to infer the inputs and interaction rules of 
their nodes from microarray data ("Kau ffmanl ri969'- 'Glass and Kauffman', 'l973l: 



[kauffman . 1993; ,Li et all . [2004: Martin et alJ . l2007 : .Samal and Jain, 2008.) . 



An n-node BN model consists of n interacting nodes which, in the context 
of biochemical networks, generally represent genes or molecular species such as 
proteins, RNAs or metabolites. Let Xi{k) denote the Boolean variable describ- 
ing the state of node i at discrete time A: = 0, 1, ... If node i represents a protein, 
then we say that if the value of the node is 1, then the protein is present, in its 
active form and its target (or substrate) present. Using De Morgan's law, the 
negation of this conjunction gives the interpretation for value 0: the protein is 
absent or inactive or its target (or substrate) absent. 

Remark 1 A more concise interpretation instead of "present" (resp. "absent") is 
"present and not being degraded" (resp. "absent or being degraded") or even "present 
and production rate greater than degradation rate" (resp. "absent or degradation rate 
greater than production rate"). 

For an n-node synchronous BN, the interactions between the nodes are mod- 
eled by a set of n Boolean interaction functions such that: 
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X,{k + 1) = F,[Xi{k), X2{k), . . . ,X„(fc)], i = 1, 2, . . . , n, 



(1) 



with Fi interaction function of node i. If the state of node i at time (fc + 1) 
depends on the state of node j at previous time fc, then node j is said to be 
an input for node i. The number of inputs of a node is called the connectivity 
of that node. The network is said to be synchronous because as expressed in 
|T|) the n nodes are updated synchronously. Also from ([T]) the dynamics of the 
network deterministic. 

Let us take an example to illustrate some essential properties of BNs. Con- 
sider the following BN: 



where symbols ^, V and A represent logical operators NOT, OR and AND 
respectively. The size of the state space of this network, i.e. the number of 
possible states, is 2" = 16 since n = A here. From ([2]), the next state for each 
possible state can be computed to obtain the state table of the network. Using 
this state table then, the state diagram can be built. This is shown in Fig[l] 
where it can be seen that the state space has been partioned into two disjoint 
sets. These sets are called basins of attraction and are denoted by 81 and B2 in 
the figure. Each basin consists of an attractor and some transient states. The 
attractor of Bl is the fixed point 1010: whatever the initial state in Bl, the 
network will converge to 1010 and stay there forever. Attractors may also be 
periodic as is the case for basin B2. The size of a basin or of an attractor is the 
number of states that constitute it. 

Now suppose that, due to random perturbations, each node of a BN has 
probability < p < 1 to switch its state (0 to 1 or vice versa) between any 
two times k and (fc + 1). Then we get a noisy BN (NBNjll. The introduction 
of disorder p is supported by the stochastic nature of intracellular processes 
coupled to the fact that, thermodynamically speaking, biochemical networks 
are open systems. 

One important difference between a BN and a NBN is that in the latter, 
transitions between the basins of the network are allowed. For example, in the 
state diagram of Fig. [l] if we perturb simultaneously nodes 2 and 3 of attractor 
state 1010 which is in Bl, then we go to transient state 1100 which belongs to 



^Here, for the sake of simplicity, perturbation probability p is supposed to be independent 
of time k and of the states of the nodes. If, for instance, for a particular node, to 1 random 
switching probability is set greater than 1 to one, then it means that random perturbations 
tend to turn the node ON rather than OFF. 



Xi(fc + 1) 

X2{k + l) 
Xz{k + l) 
Xi{k + 1) 



^Xi{k) 
^Xsik) 

Xi{k)\J X2{k)\/ Xi{k) 
Xi(fc) AX2(fc) AX3(fc), 



(2) 



B2. 
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Figure 1: State diagram of the four-node network defined by interaction rules ([2|. 
Interactions between the nodes lead to the partitioning of the state space into two 
basins Bl and B2 of respective sizes 6 and 10. 

Remark 2 Flipping the value of single nodes has been envisaged bv \Kauffmar\ \l99i ) 
in gen etic Boolean networks to study their stability to what he called minimal perturba- 
tions. \Shmulevich et al\ \200A ) proposed a model for random gene perturbations based 
on probabilistic Boolean networks (a class of Boolean networks that enclosed the class 
of synchronous Boolean networks by assigning more than one interaction function to 
each node) in which any gene of the network may flip with probability p independently of 
other genes and interpreted the perturbation events as the influence of external stimuli 
on the activity of the genes. 

Let L{x;n,p) be the probability that exactly x nodes will be perturbed 
during one time step. Then letting q = 1 — p: 

( if 0<x<n, 

L{x;n,p) = <^ g" if a; = 0, (3) 

[ if a; = n. 

This is the binomial probability distribution with parameters n and p. On 
average, np nodes are perturbed at each time step. The probability that at 
least one node will be perturbed during one time step is therefore: 

n 

^L(x;n,p) = l-(z"=r. (4) 

x = l 

Since we are mainly interested in the behaviour of the network as p varies, in the 
following n and x will be considered as parameters while p will be considered 
as a variable. Thus we should write L{p; n, x) instead of L{x; 

'^L(p\n,x) is a function of p with parameters n and x. It is not a probability density 
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p 





0.002 


0.02 


0.2 


0.5 


L{p■A^) 


0.0080 


0.0753 


0.4096 


0.2500 


Lip-A,2) 


0.0000 


0.0023 


0.1536 


0.3750 


L(p;4,3) 


0.0000 


0.0000 


0.0256 


0.2500 


L(p;4,4) 


0.0000 


0.0000 


0.0016 


0.0625 


L(p;8,l) 


0.0158 


0.1389 


0.3355 


0.0312 


L(p;8,2) 


0.0001 


0.0099 


0.2936 


0.1094 


L(p;8,3) 


0.0000 


0.0004 


0.1468 


0.2188 


L(p;8,4) 


0.0000 


0.0000 


0.0459 


0.2734 



Table 1: Some values of the function L{p; n, x). 



For < X < n, L{p;n,x) has a maximum at p = x/n and L(0;n,x) 
L{1; n, x) = 0. 

The Maclaurin series of L{p; n, x) up to order 2 is: 



L(j);n,x) 



l-np + n{n-l)p'^/2 + ... if a; = 0, 

np — n{n — l)p^ + . . . if a; = 1, 

n{n-l)p^/2 + ... if x = 2, 

+ ... if 2 <x <n. 



Thus one has: 



1 — np + o{p) if x = 0, 
L(j);n,x) ~ ^ np + o{p) if a; = 1, (5) 

o{p) if 2 < X < n. 

This means that, n being fixed, for sufficiently small p the probability that 
2 < a; < n nodes be perturbed during one time step is negligible compared to 
the probability that just one node be perturbed. Table [T] gives some values of 
n, x) rounded to four decimal places. We see that at p = 0.002, L{p; 4, 1) = 
0.0080 4p and L{p; 8, 1) = 0.0158 « Sp. 

1.2 The mean specific path 

Consider a series of Bernoulli trials where each time step defines one trial and 
where a success means that at least one node has been perturbed. From ([4]), 
the probability that the first success will occur on trial i is: 

p,^r{l-ry-\ i-1,2,..., (6) 

which is a geometric distribution with parameter r. The first moment of this 
distribution (the mean time between two successes) is: 



function. For fixed p, L(p\ n, x) is a probability. 
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r 



1 



We see that r is a decreasing function of p and n, that tends to 1 as p — > 1 and 
to oo as p ^ 0. Since time step is unity, r equals the mean number d of state 
transitions between two successes. By analogy with the mean free path of a 
particle in physic^, d will be called mean specific path. The term "specific" 
is used to recall that between two successes, the trajectory in the state space is 
specific to node interactions, i.e., it is entirely determined by node interactions. 

For sufficiently small p one has r = 1 — q" np ^ L(ji; n, 1), where the last 
approximation comes from ([5|. Hence, for sufficiently small p we may write: 



Depending on the value of p, different dynamical regimes are possible. For 
p = 0, the network is trapped by an attractor where it stays forever. For 
< p < 1, transitions between the basins of the network occur. The more 
p is close to one, the shorter the times spent on the attractors, the more the 
network suffers from functional instability. In the low p regime, the network 
may both maintain a specific activity for a long time period and change its 
activity: functional stability and fiexibility (or diversity) coexist. 

1.3 Time evolution equation 

A NBN is in fact a discrete-time Markov chain {X/j, A: = 0, 1, . . .}, where is 
the random variable representing the state of the network at time k. The state 
space of an n-node BN will be denoted by {1,2,..., E}, with E = 2" and with 
state i corresponding to binary representation of (i — 1). 

Let TTij = Pr{Xfc+i = j|Xfc = i} > be the conditional probability that the 
network will be in state j at {k + 1) while in state i at k. The matrix of size E 
whose elements are the tt^'s will be denoted by 11 and is called the transition 
probability matrix in Markov theory. The sum of elements in each row of 11 is 

(k) 

unity. Let be the probability that the network will be in state i at time k 
and denote by z^'^^ the vector whose elements are the state probabilities z^'^K 
Given an initial state p robability vector z^°\ the vectors z^^^z^^^ . . . are found 



Pi « np{l ~ np) 



i = l,2,... 




zC^+D^z^n, fc-0,1,... 



(8) 



Matrix 11 is the sum of two matrices: 



1. The perturbation matrix H', whose (i, j)th element is equal to 




if i 7^ 
if i = j, 



■'In kinetic theory of gases, the mean free path of a gas molecule is the mean distance 
traveled by the molecule between two successive collisions. 



where hij refers to («,j)th element of Hamming distance matrix H, a 
symmetric matrix that depends only on n. Element hij is equal to the 
number of bits that differ in the Boolean representations of states i and 
j. Diagonal elements of H are therefore 0. For example, if n = 2, iJ takes 
the form: 





/ 


1 


1 


M 


H = 


1 





2 


1 


1 


2 





1 




V 2 


1 


1 


/ 


of H, integ 


ers 


X — 


1,2 



(9) 



. , n appear respectively 
(") times. Thus from ([3]) and (0]), the sum of elements in each row of 11' 
must be r. 

The interaction matrix 11". If the node interactions are such that starting 
in state i the network is forced to transition to state j in one time step, 
then («, j)th element of 11" equals g" = 1 — r, otherwise it is 0. Notice 
that if ith diagonal element of 11" equals g" then state i is a fixed point 
(attractor of size 1). 



Therefore, for fixed n, any probability tt^ will be either or a function of p. 



1.4 Stationary probability distribution 



Since p does not depend on k, the TTg 's are independent o f time. The chain 
is thus homogeneous. As shown by Shmulevich et al. ( 2002f ) for genetic proba- 
bilistic Boolean networks, for < p < 1, the chain is irreducible and aperiodic. 
Consequently, for gi ven p there is a unique stationary distribution z which is 
independent of z(°) 
two equations: 



(|Kleinrocld . Il975f ) . The stationary distribution satisfies the 



z = zII and 



(10) 



In addition, we have: 



lim n'' 



n, 



(11) 



where each row of ft is equal to vector z. 

Fig. [2] shows the stationary state probabilities Zi for the network of Fig. [T] 
and two p values. Blue points correspond to p = 0.04 and red ones to p = 0.4. 
As can be seen in the figure, when p — 0.04, the stationary probabilities of the 
transient states are small compared to those of the attractor states (attractor 
states are represented in bold type in the figure). Therefore in the low p regime, 
the activity of the network in the long term is governed mostly by its attractor^. 



■*Equivalently, in the low p regime, it is mostly attractors that determine the fate of the 
cells. 
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Figure 2: Stationary state probabilities for the network shown in Fig. [T] and two 
p values. Blue points: p = 0.04; red points: p = 0.4. Attractor states are in 
bold type. The stationary probability of 81 (obtained by summing the stationary 
probabilities of the states of Bl) is 0.3270 when p = 0.04 and 0.3708 when p = 0.4. 
Since the size of Bl is 6, this probability converges to 6/16 = 0.375 as p tends to 
1. 



In fact, the stationary probabilities of the transient states tend to as p tends 
to 0. Also notice for p = 0.04 the stationary probabilities of the second attractor 
almost equal. 

As p tends to unity, the stationary probabilities tend to be uniform, i.e. as 
p ^ 1 one has zi 1/E = 1/16 = 0.0625 \fi (see the case p 0.4). Thus for 
p sufficiently close to unity, the stationary probabilities of the transient states 
are not negligible compared to those of the attractor states. 

Remark 3 For sufficiently small p, the stationary probability of an attractor state 
may be approximated by Z* /A where Z* is the limit as p ^ of the stationary occupa- 
tion probability Z of the basin containing the attractor and A the size of the attractor 
(see further 4-3)- Thus in the case of the network of Fig. QJ we get for attractor 
state 11 (fixed point) that z\\ ~ 1/3 = linip^o ^/l, and for the states of the second 
attractor, that each stationary state probability ~ 1/6 = limp_,o Z/4. 
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2 Method for the calculation of the sojourn time 
distribution in a basin of a NBN 



Consider a BN having at least two basins. Consider a basin B of size B of this 
network and a fixed perturbation probabiHty p. Let 

The square B x B sub-matrix Q contains the one-step transition probabiHties 
TTij between the states of B. For example, in the case of basin Bl of Fig. [H 
element (1,4) of Q is equal to the one-step transition probability between state 
2 (0001) and state 10 (1001). Element i of vector a represents the one-step 
transition probability between state z S B and an absorbing stat^ regrouping 
states j ^ B, i.e. 

3<^B i^B j^B 

is a row vector of length B with all elements zero and 1 is a scalar. 

From matrix Q, we can calculate the probability W^''^ that the network will 
be in basin B at time k given an initial probability vector b^"^ of length B with 
sum of elements 1: 



b(fe+i) =b(*'-)Q, fc = 0,l,..., (13) 



and by definition: 



i=l 

with = 1. 

Now we want the sojourn time distribution in B. We shall denote by S the 
discrete random variable representing the sojourn time in B, with sample space 
s = 1,2,..., and by V's the probability distribution of S, i.e. tps = Pr{S = s}. 
For a given basin, tps depends on p and initial vector b^^^. It is given by: 

i/;^ = M^(-^-i) - s = l,2,..., (14) 

while the cumulative distribution function of S is: 

i>s ^ Pr{S < s} ^ 1 ~ W^'\ s = l,2,... (15) 
Matrix Ilf, in (fT2|) represents an absorbing Markov chain ()Snelill959f ). i.e. it 



has at least one absorbing state and from every non-absorbing state (every state 
G B) one can reach an absorbing state (any state ^ B). If the chain is initially 
in state i S B, then the mean time spent in state j G B before absorption is 



'A state such that once reached, it is not possible to escape from it. 
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the (i, j)th element of fundamental matrix {I — Q) ^ I Snelll . 1959I ). Thus if ^ 
denotes the mean of Tps then: 



/i = bW(/-Q)-il. (16) 

Here 1 is a column vector of length B with each element 1. 

Also notice from the definitions of the mean of 5" and of V's that fj. can be 
expressed as follows: 

oo 

= ^ W^''^ = 1 + T^(i' + T^(2) + . . . 

k=0 

3 The problem of geometric approximation 

Let gs be a geometric distribution with parameter p^. — 1/^, i.e. Vs and gs 
have same mean. Consider the maximum deviation between the cumulative 
distribution functions (cdfs) of these two distributions, namely: 

r =max,5(''), (17) 

S>1 

with S^'^^ = — 5s| and gs the cdf of gs, i.e. gs = 1 — (1 — s — 1,2,... In 
probability theory, S* is called the Kolmogorov metric. As 2ps, S* is not defined 
for p = and, for a given basin, depends on p and b^^). It comes from ifTS]) 
that: 

6^'^ = \{l~p^y -W^'^l s = l,2,... (18) 

Thus represents the absolute error between W^^^ and its geometric approx- 
imation (1 ~ Pi_,Y . For basins of size B = 1 (one fixed point and no transient 
states), t/i^ is given by ((6]) which is geometric. Hence for such basins 5^^^ — 
Vs,p and thus 6* = Q \/p. This is a particular case and in the following we shall 
study the behaviour of S* as p tends to without any assumptions on ips- 
The probability to exit B during (fc, k + 1) is: 

Pe(fc,fc + l) = Pr{5' = fc + l|5'>fc} = l- ^^^^ , /c = 0,l,... (19) 

This probability depends on k, p and b^"^. In particular, for B = 1, pe{k,k+ 1) 
is constant and equal to r. When p > 0, any state i € B communicates with 
any state j ^ B, thus Pe{k, /c + 1) > Vfc = 0, 1, . . . and therefore probability 
W^''^ is strictly decreasing. 

Inversely, since W^*^^ = 1, we have: 

s-l 

W^'^ = Y[[1-Pe{k,k+1)], s = l,2,... (20) 

k=0 
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and thus 

s-l 

S^'^ = \{1-P^r -I[[l-Pe{k,k+1)]\, 5=1,2,... (21) 

We see that if Pe(fc, A; + 1) = Vfc then S* = 0. Also note the following: 
1. 

lim^ = cx) Vb(°). (22) 

Whatever the initial conditions, as p becomes smaller, it takes on average 
more and more time to leave the basin. Also this means thatVb(o): -> 
as p — > 0. 

2. 

limpe(fc,/e + 1) = Vfc,b("). (23) 

Whatever k and the initial conditions, the probability to leave the basin 
during (fc, k+1) goes to as p — > 0. Indeed, from equations (fT3|) and (fT9|) . 
it comes that 



Pe(fc,fc+l) = ^a.6f' =^X^r^pV-"6f , A: = 0,1,..., (24) 

ieB iGB x=l 

with rf > the number of ways of leaving B by perturbing x bits of state 



b. 



b 



(k) 



I.e. 



^6f)=l, fc = 0,l,. 



iGB 

In order to show that Vb'^"' any S^^^ tends to as p tends to 0, we introduce 
two propositions. 

Proposition 1 For any basin of a noisy Boolean network with fixed < p < 1, 
the fundamental matrix {I — Q)^^ has a real simple eigenvalue A* > 1 which is 
greater in modulus than any other eigenvalue modulus, that is X* > |A| for any 
other eigenvalue A of {I — Q)~^ ■ We have: 

lim Pe{k, k + 1) = —. 

k — ^oo A'^ 



^This number is equal to the number of elements in row i of the Hamming distance matrix 
H that are equal to x and whose column index j is such that j ^ B. 
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For a given basin, the rate of convergence depends on p and h^^^ while the limit 
1/A* depends only on p. 

From that proposition, we see that if V's is geometric then = 1/A*. For 
example, if i? = 1 then p^ ^ r ~ 1/A*. 
Demonstration. 
We have to show that: 

Um = Ah, 

with Ab = 1 — 1/A* a real simple eigenvalue of Q greater in modulus than any 
other eigenvalue modulus and < A;, < 1. 

Matrix Q is nonnegative. It is irreducible and aperiodic thus primitive. From 
the Perron-Frobenius theorem, we have that: (1) Q has a real eigenvalue A^ > 
which is greater in modulus than any other eigenvalue modulus. Since Q is 
substochastic, we have < Ab < 1. (2) A;, is a simp le root of the characteristic 
equation of Q. Moreover I Douglas and Brianl . 19991 ): 



r T 
nm — -j- = vu 

k^co A 



with V and u right and left eigenvectors associated with Ah chosen in such a way 
that M > 0, ?; > and u^v = 1. The greater k, the better the approximation 
Q*^ « X^vu^. Thus: 

Therefore: 

Xh as fc oo. 



Proposition 2 ^ is asymptotically equivalent to A*; 

lim ii = 1. 

P^o X* 

Since A* does not depend on b^*'^ the last statement is equivalent to say 
that as p tends to 0, ^ is less and less dependent on b'"' . In other words, from 
ifTB]). the elements of vector (/ — Q)~^l tends to be equal as p ^ 0. 

Demonstration. 

If ips is geometric then ^ does not depend on b'"' . Thus from (fT6|) it comes 
that: 

{I-Q)-H = til, 

or equivalently 

Ql^(l--)1. (25) 
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This means that (1 — is an eigenvalue of Q with associated eigenvector 
1. Since Q is nonnegative primitive, by the Perron-Frobenius theorem, (1 — 
must equal Ab, which implies ^ > 1. Now Afc = 1 — 1/A*. Thus ^ — \* . From 
l(25|) we see that the sum of elements in any row of Q is a constant. This is 
because when V's is geometric, the probability to leave the basin from any state 
is constant. 

When Q is not geometric, the Perron-Frobenius theorem gives: 

Qv = \bv, < Afc < 1. 

Now as p tends to 0, Q tends to a stochastic matrix (because vector a tends to 
vector null). Therefore as p tends to 0, Ah must tend to 1 and v must tend to 
1. Thus for sufficiently small p we may write: 

Ql « Afel, (26) 

or equivalently 

{I -Q)-H^\*\. 

Therefore 

M = b(°)(/ - Q)-^ « b(")A*l A*. 

Thus as p — > 0, ^ will be less and less dependent on initial conditions and the 
error in the above approximation will tend to 0. Note that since Af, ^ 1 as 
p — > we must have /i oo as p ^ which is result (|22|) . 

Remark 4 As will he discussed later, approximation /i ~ A* rrtay he good in some 
neighborhood of some p (typically in the neighborhood of p = 0.5, see Fig. 13 in 4-3)- 

From Proposition 2 now, (1 — p^J" will tend to (1 — ^/^J'Y ~ as p ^ 0. 
On the other hand, from (|26| . we get 

Qn«Agl, 

and thus: 

M^(^) =b('')l = b('"Q''l wA^ 

Hence from ifTS]) any J^"*' will tend to as p tends to 0, whatever b^''^. Thus the 
maximum will do so: 

lim(5*=0 Vb'"'. 

Since the maximum deviation between the two cdfs of i/'s and gs vanishes as 
p ^ 0, these two cdfs tend to coincide as p tends to 0. Also this means that 
(1 —p^Y is a good approximation to within ±5* of W'^'^^ at any and every s and 
that the error can be made arbitrarily small by taking p sufficiently close to 
(but not equal to 0). 
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In addition, from ([22)) we may write for sufficiently small p: 

9s = Pr{5 < s} = 1 - e'''"^!-^") « 1 - e^f"". 
Yet the cdf of an exponential distribution with parameter 7 is: 

Pr{T <t} = l- e-'^K 

For t = a and 7 = p^, the two probabilities are equal. Therefore, as p — *■ and 
Vb(°) , the cdf of V's tends to coincide with the cdf of an exponential distribution. 
Since [i is less and less dependent on b(°) as p ^ (see Proposition 2), one 
has for sufficiently small p that an exponential distribution can be used to 
approximate the time spent in a basin of a NBN. 

Remark 5 Depending on b^"-*, tps may be very close to a geometric distribution while 
p is not small. See Examples 1 and 2 below. 

Suppose that for each basin of a NBN the sojourn time S is sufficiently 
memoryless. Then, under some additional assumptions that will be discussed 
in section 5, one may define another discrete-time homogeneous Markov chain 
Yfc whose states are the basins of the network. (1) This new stochastic process 
is coarser than since only transitions between the basins of the network are 
described. (2) The size of its state space is in general much smaller than E, 
the size of the state space of the original process. (3) It has to be viewed as an 
approximation. The tilde symbol in Y^- is used to recall that in general, a basin 
does not retain the Markov property (i.e., is not geometric). 

Example 1 Calculation of sojourn time distributions and illustration of the geo- 
metric approximation problem. Let us take basin Bl of Fig. [T] with p = 0.02 and 
bf^ = 1/6 Vi = 1, 2, . . . , 6, and let us compute ipa and Qs. The two distributions are 
plotted in Fig. [3ti. The circles stand for the probabilities ills {n = 13.3530), the points 
for the probabilities Qs (p^ = = 0.0749). Fig. [Sb shows W^*'^ and its geometric 
approximation (1 — Pp)'° versus time k. Maximum deviation S* between the two is 
2.3007%. 

Table [2] gives /i and 5* for different p values and different initial conditions for the 
two basins of Fig. [T] Each of the four columns below ^ and 5* corresponds to one p 
value, from left to right: p = 0.002, 0.02, 0.2 and 0.8. The mean values of the columns 
are also given (see the rows with "Mean"). For the calculation of /x, we considered two 
types of initial conditions: (1) the network starts in one state of the basin (successively 
each state of the basin was taken as initial state) and (2) bf^ = 1/B Vi — 1,2, . . . , B 
(see "Unif." in the Table). We see that for the first three values of p the mean sojour 
times in Bl are smaller than those in B2 (see the first three columns below jj.). Thus Bl 
is less stable than B2. Also calculated for each basin and the four p values the variation 
coefficient of the sojourn times obtained with the first type of initial conditions. The 
results are (columns 1, 2, 3 and 4): 0.1876, 1.7880, 9.3454 and 15.9039% for Bl; 0.1566, 
1.3030, 2.5152 and 17.0044% for B2. Therefore ^ is less and less sensitive to the initial 
state as p decreases, which is in accordance with Proposition 2. 

Another important observation in this table is that, depending on the initial condi- 
tions, V's may be close to a geometric distribution while p is not small (see in the table 
the values of 5* when p = 0.2). Fig. |4] shows the functions S*(p) for the two basins of 
Fig. [Hand the uniform initial condition. We see that in both cases <S*(p) — > when 
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p ^ 0, 5* (0.5) = {tps is geometric at p = O.sj^ and 5*{p) tends to a global maximum 
as p — > 1. For basin Bl (blue curve), there is one local maximum at p = 0.112, while 
for basin B2 (green curve) there are two local maxima, one at p = 0.0530 and the 
other at p = 0.3410, and one local minimum at p = 0.2600 (0.0677%). As mentioned 
above, function 5*(p) depends on initial conditions. For example, if the network starts 
in state 0000, then 82 has still one local minimum but this now occurs at p = 0.3640 
(0.0830%). 

Example 2 In this example, we illustrate the fact that while p is not small, ips may 
be close to a geometric distribution, depending on b'-''^ We take basin Bl of Fig. [T] 
and calculate tps when initially the network is (a) in state 0101 and (b) in state 0110. 
Fig. [5]shows the results. As can be seen in the figure, when the network starts in 0110, 
ips is close to a geometric distribution which is not the case when the initial state is 
0101. 

To compare the two approximations, we computed the variances and of ips 
and (?s respectively and the total variation distance drv between ips and gs (another 
probability metric) which is given by: 



dTV = ^'^li's ~ 9s 



2 

The values of these parameters are given in Table [3] We see that the total variation 
distance is about 10 times greater when the initial state is 0101 than when it is 0110. 
The more geometric ips , the smaller drv , the better the approximations ag and 
fi « A*. 

Till now, we have assumed each node of the network has the same probabil- 
ity of being perturbed. In Example 3 below, we look at how the network of Fig. 
[U behaves when one of its nodes is perturbed with a probability which is high 
compared to the other nodes. Suppose nodes represent proteins. Within the 
cell interior, some proteins may be more subject to conipeting reactions than 
others. These reactions may be assumed to act randoml}!^, either negatively or 
positively, on the state of the target proteins. For example, some reactions may 
lead to protein unfolding while oth ers, Hke t hose involving molecular chaper- 
ones, may rescue unfolded proteins (Dobson, 2003h . On the other hand, some 



proteins may be more sensitive to physico-chemical factors, like temperature or 
pH, increasing their probability to be perturbed. 

Example 3 In this example, it is assumed that one node is perturbed with a proba- 
bility which is high compared to the other nodes. The behaviour of the network when 
p — > 1 is complex and will not be discussed here so we take pi = P2 = Pa = 10^'' 
(the first three nodes are rarely perturbed) and < p4 < 0.9. Fig. [6] shows z for two 
values of p4. For p4 = 0.04 (the blue points) the network behaves as if all the nodes 
had the same probability p = 0.04 of being perturbed (see Fig. [2] the blue points). 
Simulations have shown this to be true whatever the three rarely perturbed nodes. 
When p4 increases in ]0, 0.9], the stationary state probabilities do not tend to be equal 
(see Fig. [6l the case P4 = 0.9), rather, some transient states, typically state 0011 of 
82, tend to be more populated at the expense of attractor states (see attractor states 



''Whatever the initial conditions, at p = 0.5, n = 8/5 for Bl and 8/3 for B2. See 5.1.2 for 
analytical expression of fi when p = 0.5. 

*Due to the fluctuating nature of intracellular processes. 
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Figure 3: Sojourn time distribution and geometric approximation for basin Bl 
of Fig. [1] with p = 0.02 and = 1/6 Vi = 1,2,..., 6. (a) Circles: sojourn 
time distribution ijjs (with mean /i = 13.3530); points: approximating geometric 
distribution gs (with parameter = l//i = 0.0749). (b) Circles: W^'^'^; points: 
geometric approximation (1 — p^j)'^. Maximum deviation 5* is 2.3007%. 
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Figure 4: Maximum deviation 6* (in %) as a function of p for the two basins of 
Fig. [T] and initial conditions 6-"' = 1/_B Vi = 1, 2, . . . , S. Blue curve: basin Bl. 
Green curve: basin B2. One has 5* {0.5) = for both basins (-ips is geometric at 
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Ini. cond. /i S* 



Bl 


0001 


126.00 


13.52 


2.37 


1.54 


0.39 


3.44 


7.47 


10.39 




0101 


126.00 


13.52 


2.39 


2.07 


0.39 


3.47 


8.49 


4.09 




0110 


125.50 


13.02 


1.94 


2.28 


0.00 


0.00 


0.92 


10.66 




1001 


126.00 


13.52 


2.39 


2.07 


0.39 


3.47 


8.49 


4.09 




1010 


125.50 


13.02 


1.94 


2.28 


0.00 


0.00 


0.92 


10.66 




1101 


126.00 


13.52 


2.37 


1.54 


0.39 


3.44 


7.47 


10.39 




Mean 


125.83 


13.35 


2.23 


1.96 


0.26 


2.30 


5.63 


8.38 




Unif. 


125.83 


13.35 


2.23 


1.96 


0.26 


2.30 


4.59 


2.72 


B2 


0000 


252.85 


27.68 


4.27 


3.77 


0.39 


2.96 


2.31 


5.40 




0010 


251.86 


26.81 


4.08 


2.25 


0.00 


0.11 


0.50 


20.87 




0011 


252.36 


27.26 


4.22 


3.77 


0.20 


1.63 


2.59 


5.42 




0100 


251.86 


26.80 


4.02 


3.60 


0.00 


0.12 


1.99 


0.93 




0111 


251.86 


26.80 


4.02 


3.60 


0.00 


0.12 


1.99 


0.93 




1000 


251.86 


26.80 


4.02 


3.60 


0.00 


0.12 


1.99 


0.93 




1011 


251.86 


26.80 


4.02 


3.60 


0.00 


0.12 


1.99 


0.93 




1100 


252.36 


27.26 


4.22 


3.77 


0.20 


1.63 


2.59 


5.42 




1110 


251.86 


26.81 


4.08 


2.25 


0.00 


0.11 


0.50 


20.87 




nil 


252.85 


27.68 


4.27 


3.77 


0.39 


2.96 


2.31 


5.40 




Mean 


252.16 


27.07 


4.12 


3.40 


0.12 


0.99 


1.88 


6.71 




Unif. 


252.16 


27.07 


4.12 


3.40 


0.12 


0.81 


0.24 


2.81 



Table 2: Mean sojourn time /i and maximum deviation S* (in %) for the two basins 
of Fig. [H Two types of initial conditions have been considered: (1) the network 
starts in one state of the basin (successively each state of the basin is taken as initial 
state) and (2) each element of b'^*'' is one divided by the size of the basin (see "Unif." 
in the Table). The four columns below a parameter (/i or S*) correspond from left 
to right to p = 0.002, 0.02, 0.2 and 0.8. The mean value for each column is also 
indicated (see "Mean"). Attractor states are in bold type. 





0101 


0110 


S* (%) 


8.49 


0.92 


dTV (%) 


10.75 


1.18 




2.38 


1.91 


-I 


3.34 


1.82 


X* 


2.00 


2.00 




2.39 


1.94 



Table 3: Comparison between two geometric approximations. See Fig.[5]for details. 
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Figure 5: Comparison between two geometric approximations. Two different initial 
states are taken in basin Bl of Fig. [T]with p ~ 0.2 (see table[2|. Circles: probabil- 
ities ips- Points: geometric approximation gg. (a) Initial state 0101. S* = 8.49%, 
fi = 2.39. (b) Initial state 0110. 6* = 0.92%, n = 1.94. 
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Boolean states 



Figure 6: Stationary state probabilities for the network shown in Fig. [T]and non- 
equiprobable node perturbations. The first three nodes are rarely perturbed: pi = 
P2 = P3 = 10^'^. Blue points: p4 ~ 0.04; red points: p4 — 0.9. Attractor states 
are in bold type. 



1000, 1010, 1011 and 1110 in Fig. [6lP|. 

Fig. [TK shows distribution tpa with approximating distribution gs for basin B2, 
P4 = 0.1 and state 0011 as initial state. Maximum deviation S* is 8.2011%. Note the 
splitting of ips into two subdistributions, one for odd frequencies and the other for 
even frequencies. Probabilities W'*^^ and approximations (1 — P/j)*" are plotted in Fig. 



4 The geometric approximation in (n, C) networks 

In Example 1, we illustrated, using the network of Fig. [TJ the fact that as p 
tends to 0, the cdf of the time spent in a basin of attraction approaches the cdf 
of a geometric (resp. exponential) distribution and that the expected sojourn 
time tends to be independent of initial conditions (see the decreases of S* and 

^At the cell population level, this means that a cell population would exhibit phenotypes 
that could not have been observed (or not easily observed) in the low p regime. Note that 
these phenotypes are not necessarily survivable for the cells. 
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Figure 7: Sojourn time distribution and geometric approximation for basin B2 of 
Fig.ffl The first three nodes are rarely perturbed: pi = P2 = Pz = lO^'^; pi = 0.1. 
Initially, the network is in state 0011. (a) Circles: distribution (with mean 
/Li = 22.9653); points: approximating geometric distribution Qs (with parameter 
p^ = 0.0435). For the sake of clarity, sojourn times s > 70 have been omitted, (b) 
Circles: M^'^'^'; points: geometric approximation (1 —p^)'^. Maximum deviation 5* 
is 8.2011%. 
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variation coefficient as p becomes smaller, respectively in Table [2] and in the 
text). As will be seen later, for this network, a two-state discrete-time (resp. 
continuous-time) homogeneous Markov chain may be used for approximating 
basin transitions in the low p regime, i.e., a coarse-grained description of this 
network exists in the low p regime. 

Now we address the problem of geometric approximation in randomly con- 
structed (n, C) networks. A random {n, C) network is bui lt by randomly c hoos- 
ing for each node C inputs and one interaction function l|Kaufrmanl . ll993ll . 



4.1 {n,C) networks and confidence intervals 

Six (n, C) ensembles were examined taking n = 6, 8 or 10 (i? = 64, 256 or 1024) 
and C = 2 or 5. Each basin of a randomly constructed network was perturbed 

with p = 0.002, 0.01, 0.02, 0.1, 0.3, 0.5 and 0.8 . 

According to the classification o f lKauffmanl \im± . ( n, 2) networks are com- 



plex while (n, 5) ones are chaotic. One difference between these two ensembles 
of networks, which is of particular interest here, is that in the former "// the 
stability of each state cycle attractor is probed by transient reversing of the ac- 
tivity of each element in each state of the state cycle, then for about 80 to 90 
percent of all such perturbations, the system flows back to the same state cycle. 
Thus state cycles are inherently stable to most minimal transient perturbations." 
I Kauffmanl . Il993l . p. 201). In chaotic n etworks, the stab ility of attractors to 



minimal perturbations is at best modest l|Kauffmanl . [l99i p. 198) 



For each (n, C) ensemble, we generated about 2500 basins (which corre- 
sponds to about 700 to 800 networks depending on the ensemble). We es- 
tablished confidence intervals for three statistical variables of which two are 
probabilities: 

1. Consider an n-node network. For each of its basins, one can define two 
conditional probabilities: (1) the conditional probability a of leaving the 
basin given that one attractor bit out of nA has been flipped, with A the 
size of the attractor, and (2) the conditional probability /3 of leaving the 
basin given that one basin bit out of nB has been flipped. 

In each basin sample, a small proportion of basins having a = were 
found. The 25th percentile of the relative size B/E of a = basins was 
more than 0.8 for C = 2 networks and 0.9 for C = 5 ones. Thus a = 
basins are most often big basins. We calculated 2 ratios: ratio R, between 
the median mean sojourn times of a = and a > basins and ratio 
K* between the maximum mean sojourn times of the two types of basins. 
These ratios are given in Table [4] for (n, 2) networks and four p values. In 
the case of a = basins, the computation of ips for p = 0.01 and stopping 
criterion 'ips > 0.9999 varies between a few minutes to several days using a 
PowerEdge 2950 server with 2 Quad-Core Intel Xeon processors running at 
3.0 GHz. As can be seen from this table, with p decreasing, the maximum 
mean sojourn time of a = basins increases drastically compared to that 
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p 







0.03 


0.05 


0.1 


0.5 


K n = 


6 


30.1612 


19.9012 


12.1390 


6.8534 


n = 


8 


44.3541 


26.0081 


15.1438 


11.9929 


n ~ 


10 


35.5351 


21.6376 


14.5090 


10.3930 


K* n = 


6 


122.3414 


19.3446 


3.6270 


1.0000 


n = 


8 


300.3482 


32.4353 


3.1643 


1.0000 


n ~ 


10 


542.8109 


67.3100 


4.9887 


1.0000 



Table 4: Comparison between median mean sojourn times of a = and a > 
basins (ratio k) and between maximum mean sojourn times of both types of basins 
(ratio K*) for {n,2) networks and four p values. Calculation of ^ps- at time 0, the 
states of the basin are equiprobable. The proportions of a = basins for samples 
n = 6, 8 and 10 were found to be 4.0759, 3.1989 and 2.5971% respectively. 



of a > basing. The proportion of a = basins varies in our samples 
from 1.7 to 4.1% depending on the ensemble. For fixed C, it is a decreasing 
function of n and for fixed n, it is smaller in the chaotic regime than in the 
complex one. Although we do believe that a = basins have not to be 
rejected from a biological interpretation perspective, networks with at least 
one a = basin were omitted during the sampling procedure (essentially 
because of the high computational time that is needed to compute i/js 
when a = and p is small) . 

Fig. [8] shows the conditional probabilities a for sample (8,2). For an 
n-node network, one can define n main probability levels 1, 2, . . . , n corre- 
sponding to probabilities 1/n, 2/n, . . . , 1. The percentage of data points 
located in main probability levels are given in Table [5] for the six samples 
and the two conditional probabilities a and /3. It can be seen from this 
table that whatever the connectivity and the conditional probability, this 
percentage decreases when n increases. For a given conditional probability 
and any n, it is greater in the complex regime than in the chaotic one. 

The statistical analysis of the six data samples has shown the following. 
For C = 2, there is a clear quantization of a and a weak one of f3. Also 
the histograms of a and (3 are symmetric (skewness equal to ^ 0.2). For 
C = 5, the quantization of a is weak and there is no obvious quantization 
of p. The histograms of a and f3 are negatively skewed (skewness equal to 
~ —0.9) and the last two main probability levels (n — l) and n are strongly 
populated compared to the other ones. 

^"For = basins, the conditional probability of leaving the basin given that 2 < x < n 
bits of an attractor state have been perturbed may be non null. However L{p;n,x) = o{p) for 

2 < X < n. On the other hand, for transient states and sufficiently large k one has Sf ^ 
as p ^ 0. 
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Figure 8: Conditional probabilities a > sorted in ascending order, sample (8,2). 
About 74% of the data points are in main probability levels. These correspond to 
probabilities 1/8,2/8,...,!. 



c 


n 


a 


/3 


2 


6 


78.6 


52.0 




8 


73.9 


45.9 




10 


70.2 


35.3 


5 


6 


56.6 


31.5 




8 


50.8 


23.5 




10 


43.6 


15.4 



Table 5: Percentages of data points located in main probability levels for the six 
ensembles {n,C) and the two conditional probabilities a > and (3 > 0. 



24 









a 


/3 




c 


n 


Mean 


Median 


Mean 


Median 


2 


6 


0.5057 


0.5000 


0.4817 


0.5000 




8 


0.4523 


0.4375 


0.4294 


0.4144 




10 


0.4149 


0.4000 


0.3907 


0.4000 


5 


6 


0.6869 


0.7917 


0.6852 


0.7778 




8 


0.7032 


0.8125 


0.7038 


0.8219 




10 


0.7050 


0.8125 


0.7106 


0.8375 



Table 6: Mean and median of conditional probabilities a > and /? > for the six 
ensembles {n, C). 



Table [6] gives the mean and median of probabilities a and [i for the six 
ensembles (n, C). For fixed n, the mean and median of both conditional 
probabilities are higher for chaotic basins than for complex ones. Addi- 
tionally, for C = 2 basins, the mean and median decrease with n while for 
C = 5, both increase with n (except for the median of a). 

2. The third statistical variable is the mean time spent in a basin, /i. For 
C ~2 and any n, there is a clear quantization at p = 0.002 which rapidly 
disappears as p increases. For C = 5 and any n, there are only two levels 
of quantization at p = 0.002 and these rapidly vanish as p increases. Also 
notice that the minimum of ^ is equal to the mean specific path d = t 
which, according to does not depend on C. 

Most of the limits of the 95% confidence intervals ranged between 1.5 and 
2.5% (for a and /3: confidence interval for the mean if C = 2, for the median 
and for the trimmed mean if C = 5; for /j,: confidence interval for the trimmed 
meani"!. 

4.2 Simulation results and discussion 

Most of the basin variables (such as /i or (5*) were positively skewed. Therefore 
for each of these variables we calculated quartiles Qi, Q2 (the median) and Qa- 
For the calculation of ^, two types of initial conditions were considered: 

1. Each element of b^^^ is equal to 1/B. We call this condition the uniform 
initial condition. 

^^Two methods were used: a nonparametric method based on the binomial distribution (for 
the median only) and the bootstrap method (used for the median and the trimmed mean). 
For the trimmed mean, we averaged the sample data that were (1) between the 25th and 75th 
percentiles, (2) less than the 75th percentile and (3) less than the 90th percentile. 
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2. We pick one state of the basin at random and place initially the network 
in that state (the elements of b'^"^ are all except one which is equal to 
1). We call this condition the random initial condition. 

Let's start with the basin variable /i. First the uniform initial condition. 
CoefHcients of skewness for ^ distributions mostly ranged from 5 to 15 with 
mean of 10.0568 (strong skewness). For example, for sample (8, 5) with p = 0.01, 
we obtained a coefRcient of skewness of 10.2415, a mean of 44.7471, a median 
of 15.7922, a 75th percentile of 26.0523, a 99th percentile of 608.4743 and a 
maximum of 3.0117 x 10^. A small proportion of the mean sojourn times are 
therefore very far from the median (taken here as the central tendency). The 
three quartiles of ^ versus Inp are shown in Fig. [9]for the six ensembles (n, C). 
The blue squares correspond to C = 2 and the red triangles to C = 5. The 
size of a symbol is proportional to ?i. For the sake of clarity, the abscissae of 
the points corresponding to ensembles (6, C) and (10, C) have been translated 
(respectively to the left and to the right of the p values) . 

It can be seen from Fig. [9] that: (1) for fixed p and any C, the median of ^ 
is a decreasing function of the size n of the network (not true at p = 0.5 and 
0.8), although a decreases with n (see Table [6]). (2) For fixed C and any n, the 
median of /i is a decreasing function of p. (3) For fixed p and any n, median /i of 
chaotic basins (C = 5) is less than that of complex ones (C = 2). At p = 0.002, 
the medians of ^ for C = 5 ensembles are respectively 1.58, 1.86 and 2.02 times 
smaller than those for C = 2 ones. Similar but smaller v alues were f ound for 
p = 0.01 and p = 0.02. We therefore confirm the results of Kauffman ( 19931 . p. 



198, 201 and 488-491) that chaotic basins are less stable to node perturbations 
than complex ones. 

We looked for the relationship between the median of fj. and p. We found 
that for sufficiently small p: 

Qi2;,) « |, (27) 

i.e. the median of the mean sojourn time is inversely proportional to p. The 
proportionality constant C2 has been estimated by the least squares method 
for the six ensembles (n, C) and was found to decrease when n increases (only 
the first three data points were fitted, i.e. the points with abscissa p — 0.002, 
0.01 and 0.02). For C = 2, we obtained C2 = 0.3342, 0.2871 and 0.2510; for 
C = 5, C2 = 0.2117, 0.1547 and 0.1242. Thus for fixed n, chaotic basins are less 
sensitive to a variation in p than complex ones. The result of the least squares 
fit is presented in Fig. [lOl Graph (a) shows the hyperboHc relationship between 
the median of and p for the six ensembles (n, C). For the sake of clarity, the 
functions were drawn up to p = 0.05. When a logarithmic scale for each axis 
is used, one obtains the graph (b) which shows for the three ensembles (n, 2) 
a linear relationship between lnQ(2;At) and lap at low p (up to p = 0.02 on 
the graphjll- The same rule applies to the three ensembles (n, 5). We will see 
further how to express C2 in function of n and median a probability. 



^The three straight lines in Fig. llOb were obtained by the least squares method applied to 
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Figure 9: Median and interquartile range of mean sojourn time /^i versus Inp for the 
six ensembles (n, C). Calculation of tps with the uniform initial condition. Quartiles 
Qi and Qs are indicated by horizontal bars. Blue squares: complex regime (C = 2); 
red triangles: chaotic regime (C = 5). The size of a symbol (square or triangle) is 
proportional to n. 
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Figure 10: Fit of ([27]) by the least squares technique, (a) Linear scale for both axis, 
fits for the six ensembles (n, C). For the sake of clarity, the points with p > 0.05 
have been omitted. Blue squares: C = 2; red triangles: C = 5. The size of a 
symbol is proportional to n. (b) Logarithmic scale for both axis (In-ln), fits for 
ensembles (n, 2) only. The size of a symbol is proportional to n. 
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With the random initial condition, the quartiles of fj, are very close to those 
obtained with the uniform initial condition. 

We then calculated for each basin the relative error between the mean so- 
journ time obtained from the uniform initial condition and that obtained from 
the random initial condition. This error is null at p = 0.5 (as for the network 
of Fig. [U this is because ips is geometric at p = 0.5). The quartiles of the error 
tend to as p — > and they reach a maximum at p = 0.1 whatever the network 
ensemble. Thus the smaller p, the more /i is independent of initial conditions. 

To end, we turn to the maximum deviation S* . First the uniform initial 
condition. The results for the six (n, C) ensembles are presented in Fig. [iTJ At 
p = 0.002 and for fixed C, the quartiles increase linearly with n, except the first 
quartile of C = 2 basins which is constant and approximately equal to 0.0005% 
(whatever n, 25% of C = 2 basins have a sojourn time that is geometric or 
closely follows a geometric distribution). 

For both connectivities, the functions Q(2;S*){p) and Q(3;5*)(p) have similar- 
ities with the functions S*{p) of Fig. |4l they tend to when p — > 0, they have 
at least one local maximum in < p < 0.5 and are null at p = 0.5. 

With the random initial condition, the second and third quartiles of S* in- 
crease for most of the six ensembles {n, C) compared to the uniform case. Qual- 
itatively, the behaviour of the quartiles with respect to p is similar to the one 
found with the uniform initial condition. 

To summarize sections 3 and 4, we have the following proposition: 

Proposition 3 Consider a basin of a noisy Boolean network. As p tends to 
then: 

1. whatever the initial conditions, /i ^ oo. 

2. /i tends to he independent of the initial conditions: fi ^ X* with A* > 1 
the Perron- Frohenius eigenvalue of fundamental matrix {I — Q)~^ . 

3. whatever the initial conditions, 5* — > 0, where 6* is the Kolmogorov dis- 
tance between the cumulative distribution function of the sojourn time S 
in the basin and that of a geometric distribution gs having the same mean 
as S. 

4- gs converges to an exponential distribution. 

Proposition 3 does not guarantee the existence of a coarser representation 
of a NBN in the low p regime. What can be said from this proposition is that 
if such a representation exists, then it is in general an approximation of the 
original process and it can always be expressed in a discrete or continuous time 
framework. Thus we are lead to the proposition below that will be discussed in 
more details in the next section: 

the linearized problem: Y2 = a2 — X where I2 = liiQ(2;/j), 0,2 = lnc2 and X = Inp. Only the 
first three data points were fitted. 
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Figure 11: Geometric approximation for (n, C) ensembles. Median and interquartile 
range of S* (in %) versus Inp. Calculation of i/'s with the uniform initial condition. 
Blue squares: C = 2; red triangles: C = 5. As p — *■ 0, the three quartiles of 6* 
tend to 0. 
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Proposition 4 Consider a noisy Boolean network and suppose there exists a 
coarser time-homogeneous representation of this network in the low p regime. 
Then the dynamics between the basins of the network may be approximated by 
a system of linear ordinary differential equations with size being equal to the 
number of attractors of the network. 

Recall that Proposition 3 does not mean that i/'s can never be geometric 
nor be approximated by a geometric distribution when p is not sufficiently close 
to 0. For example we know tps is geometric when p = 0.5 or when B = 1 
whatever < p < 1. One difference between a strongly and a weakly perturbed 
network is that in the latter, since the mean specific path d is large compared to 
1, the network spends long periods on the attractors without being perturbed. 
Under normal conditions, biochemical networks are supposed to work in the low 
p regime because as explained in 1.2, this regime is associated with functional 
stability. 

4.3 Approximation formula for the mean time spent in a 
basin of a NBN 

For sufficiently small p, (i, j)th element of 11' can be approximated by (i ^ j): 



p if hij = 1, 
if hij > 1. 



From l(24|) then we may write: 

p,{k,k + l)^pY,^}b\'\ (28) 

ieB 

We see that in the low p regime, the probability Pe{k, fc + 1) depends on time 
and initial conditions. 

Now for sufficiently small p, we have the following: 

pRil/npa, (29) 
which, from Proposition 2, is equivalent to: 

lim npaX* = 1. 

To show approximation formula l[29|) . we consider two cases: 

1. Suppose for any given < p < 1, the elements of vector a are equal. Then 
from ((24l) Peik, k+ 1) is constant which means ips geometric. Thus taking 
p sufficiently small, T} must be the same for all state i G B and therefore 
from (EHl): 



l//i = p,, « = npa, (30) 
with a = AT^ /nA and A the size of the attractor. 
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Note that: (1) if A = B = 1, it comes that = r k, np and therefore 
^ = r « 1/np. (2) \{ A ~ B = 2, then tps is geometric since in this case 
Q is a symmetric matrix of size 2. If one starts with probabiUty 1 in one 
of the two states, then each probabihty bi will be periodic with period 
2. (3) If lAs is geometric then the parameter of the distribution may be 
written: 

n 
x = l 

with ai = a and a2, as, . . . , conditional probabilities defined as a 
except that rather than perturbing one bit we perturb 2,3,...,n bits 
simultaneously. Taking p sufficiently small in this formula, we retrieve 
approximation ((30l) . 

2. If ips is not geometric, then for sufficiently small p it comes from l(28l) that 
Vb^'^^ Pe{k,k+ 1) must tend to pX^ieA^i/^ as fc oo since 6^'^'' must 
tend to a value which is close to Vi S B \ A and to a value which is close 
to 1/Ayi £ A. From Propositions 1 and 2 then, we get (|29| with: 

a^Y.^l/nA. 

Remark 6 We see from \29fl that the quantization of /i observed at p — 0.002 for 
(n, 2) ensembles (and to a lesser extent for chaotic ensembles) is a direct consequence 
of the quantization of a for these ensembles (see Fig. \^ and point 2 of 4-1)- 

From ^ and ([29|) we get for sufficiently small p that: 

^l « r/a, (31) 

Therefore, in the low p regime, the mean time spent in a basin of a NBN is 
approximately proportional to the mean specific path d — t. As shown in 
Fig. [12] for the six ensembles (n, C) (In- In plot), the relative error done in 
approximation l(3T|) decreases as p becomes smaller. It can also be seen in the 
figure that in the low p regime, C being fixed, median increases with n, 
and that while at fixed n the median error for chaotic basins is smaller than 
for complex ones, the interquartile range for the former is greater than for the 
latter. 

The statistical basin variables in equation l(29|) are fi and a. Since the hy- 
perbolic function is strictly monotone decreasing, the median of 1 /a is equal to 
the inverse of the median of a. Thus the constant C2 in equation (|27|) must be 
approximately equal to C2 = i/nQ(2-a}, where Q(2;a) stands for the median of 
a. For C = 2, the values of £2 were found to be (n = 6, 8 and 10) 0.3333, 0.2857 
and 0.2500; for C = 5, we found £2 = 0.2105, 0.1538 and 0.1231 respectively. 
These values of C2 are indeed very close to the values of C2 that have been ob- 
tained by the least squares method in subsection 4.2 (the relative error between 
C2 and £2 ranges from 0.25 to 0.90%). 

The first and third quartiles of ^ satisfy a relation of the same type as ((27|) : 
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Figure 12: Approximation of /i using ([31]) . Quartiles of relative error fr (in %) 
versus p (In-ln frame) for the six ensembles {n,C). The initial condition for the 
calculation of is the uniform one. Blue squares: C = 2; red triangles: C = 5. 
The size of a symbol (square or triangle) is proportional to n. 
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Figure 13: Comparison between /i (blue curve), A* (green curve) and 1/nap (red 
curve) versus p (In-ln plot) for basin Bl of Fig. [l]and uniform initial condition. 



where ci ^ ci — l/ri(3(3;a) and C3 ~ £3 — 1/?^(5(i;q). The first (resp. third) 
quartile of /i is thus Hnked to the third (resp. first) quartile of a. For each 
ensemble {n,C), we can define as many constants Ci as there are percentiles. 

Fig. [13] shows fi (in blue), A* (in green) and 1/nap (in red) versus p (In-ln 
plot) for basin Bl of Fig. [1] and the uniform initial condition. For small p 
{p < 10^^ on the figure), these three functions behave identically. For p = 0.5, 
we see that /i = A*. For basin B2 of Fig. [l]and the uniform initial condition, ^ 
is very close to A* whatever < p < 1 (figure not shown) . Thus formula « A* 
may be accurate even if p is not small, which is not the case for approximation 
/z « 1/nap. 

Since p^ = 1/^, from ^ and ((29l) we can express p^ through some approx- 
imation formulas valid in the low p regime: 

p^, w npa. (32) 
Since for sufficiently small p, r np, we may write: 
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Pf_i ~ ra. 

If a is in the ith main probability level (see 4.1), then: 



(33) 



Pfj, ^ ip, i = l,2,...,n. (34) 

If ips is geometric, then its parameter is equal to p^. Only in this case is the 
probability of leaving the basin during one time step equal to the inverse of the 
mean sojourn time fi. If lAs is close to a geometric distribution (for some or all 
b(°)), then Pe{k, fc + 1) is approximately constant, i.e. Peik, k + 1) p^ = 

Let us mention another property. Let Peik, k+A) = 'J2f=i Peik+j — l, k+j). 
Then 

lim pe{k, k + A)^ —. 

k—^oc A 

Equivalently, at fixed p, the mean probability of leaving a basin calculated 
over a period A of the attractor tends to be constant as time increases. The 
convergence is much more rapid than the one ofpe{k,k + l) (see Proposition 1). 

Finally, we shall establish the general expression of the mean sojourn time 
whenp = 0.5. The number of ways any state of an n-node NBN can be perturbed 
is: 

E(:)-2"-1, (35) 

where the last equality follows from the Binomial Theorem. This means that 
whatever the perturbed state, any of the (2" — 1) other states is reachable from 
that state by applying the appropriate perturbation combination out of the 
(2" — 1) possible perturbation combinations (the chain is irreducible). Hence: 

n 

^ = 2" - B Vi e B, 

x=l 

with, as before, B the size of basin B. For p = 0.5, we know that ips is geometric, 
i.e. the probability does not depend on k nor on the initial conditions. Thus: 



2" - S 

= (36) 



and thus 



. = A* = ^. (37) 

The bigger the basin, the smaller p^, the higher ^. The Perron- Frobenius eigen- 
value of Q is B/E and that of (/ - Q)-^ is E/{E - B). When B = 1 we get 
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^ ~ T = \* = E/{E ~ 1). Notice that, since formula ((37l) depends only on E 
and B, it is also valid when a = 0. 

Another way to get (|36|) is to notice that when p — 0.5, each element of 11' 
(and n") is equal to (0.5)" = l/E. This means that 

E-B 
flj = — - — Ml G B. 

hj 

Since does not depend on time nor on initial state i, ips must be geometric. 

5 Method for the reduction of a NBN 
5.1 Discrete-time reduction 

5.1.1 The low p regime 

Consider an i?-basin NBN with state space size E and Markov representation 
and suppose that the network has no a = basin. We want to find a discrete- 
time homogeneous Markov chain {Yfc, /c = 0, 1, 2, . . .} with state i of the chain 
representing basin i of the network, i.e. the state space of Yfc is {1, 2, 3, . . . , i?}. 
In the theory of Markov processes, would be called a reduced chain or 
aggregated chain because R < E. The problem of reducing a Markov chain to a 
chain with a smaller state spac e , the s o-called "state spa c e exp losion problem", 
is not new l|Kemenv and Snelll . Il960l : iFredkin and Ricd . Il986h and is still an 



active field of research in the theory of Markov proc e sses (iGuedon et al . 20061 : 



Grone et all l2008l : IWeinan et all . 120081 : Izhao et al.1 . l2009f ). Here, we do not 



need to aggregate X^. The aggregation is fixed by the interactions that occur 
between the components of the network. 

The expressions for the transition probabilities ftij = Pr{Yfc+i = j|Yfc = i} 
between the basins of the network are found as follows. For convenience, the 
validity of these formulas are discussed further below. Prom l(32|) . we write: 



TTji = 1 - npai, 

where ai denotes the a probability of basin i [ai > 0, Vi). Thus we have: 

TTij = npaj, i = 1, 2, . . . , i?. (38) 

Probability ai can be expressed as a sum of probabilites: 

Q!j = ^ ftij , (39) 

with aij the conditional probability for a transition between basins i and j to 
occur given that one bit of attractor i has been perturbed {an = 0). If aij > 
(z 7^ j), then transition probability fcij is taken to be: 

TTjj = npaij, ij^j. (40) 
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Therefore in the low p regime, Tfy is the product of two probabiHties: the proba- 
bihty that one node be perturbed during one time step, which is approximately 
equal to np for small p, and the conditional probability aij . 

Now if Uij = 0, one may transition to basin j by perturbing at least two bits 
of an attractor state or at least one bit of a transient state. Since in the low 
p regime the network is rarely found in transient states and L(ji;n,x) = o{p) 
for 2 < X < n, when aij = 0, transitions i ^ j are rare events compared to 
transitions i ^ j for which aij > 0. Therefore we take itij = 0. 

By supposing that the time spent in any basin B of a NBN is geometric with 
mean 1/npa, we neglect transitions of order 1 from transient states as well as 
transitions of order 2 or more (transitions from transient or attractor states due 
to perturbations affecting two or more nodes simultaneously). This means that 
while the original chain is irreducible, the reduced chain may not be irreducible 
anymore. If this is the case, the reduction method may give inaccurate results. 
Suppose that reduction of a NBN gives two sets of basins, each containing two 
basins that communicate with each other (each basin is accessible from the 
other) , and that those sets are closed (by perturbing any node of any attractor 
state of any set, the other set cannot be reached). In Markov theory, such sets 
are called closed communicating classes. If we start in one set with probability 
one, then the state probability in the other set calculated from the reduced 
matrix will be at any time. Now if we solve the original chain, this will not be 
the case. If the stationary probability for the initially empty set is not negligible, 
then it will take a long time to approach this probability with good accuracy 
but it will. Another difference between these two chains is that the reduced 
one has an infinity of stationary distributions while the original one a unique 
stationary distribution. 

If, starting in any basin, one can reach any other basin by applying single 
node perturbations to attractor states, then the reduced chain is irreducible. 
In this case, the smaller p, the more accurate the reduction method. If the 
reduced chain is not irreducible, the reduction method is not guaranteed to 
work properly. 

Remark 7 We investigated the case of reducible chains. Let ipij be the probability 
distribution of the time spent in basin i given bo and arrival basin j . When aij = 0, 
the first moment fiij of ipij rnay strongly depend on bo. We found some cases (some 
basins with some bo ) in the low p regime for which ^ij was significantly smaller than 

To illustrate the chain reduction method, we chose a randomly generated 
(8, 2) network having i? = 4 basins of size 72, 120, 36 and 28. The size of the 
corresponding attractors were 6, 6, 1 and 3. We considered twop values, namely 
0.01 and 0.1. The reduced matrix 11 when p = 0.01 was found to be: 

/ 0.9633 0.0200 

~ _ 0.0133 0.9700 

0.0400 

\ 0.0333 0.0200 



0.0167 \ 

0.0100 0.0067 
0.9600 

0.9467 I 



(41) 
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This matrix is irreducible, i.e. any basin is accessible from any other basin 
by applying single-node perturbations to attractor stateJ^. Also note that the 
second basin, which is the biggest one (120 states), is the only state of which 
is reachable from any other state. 

The four basin occupation probabilities versus time are shown in Fig. [HI 
The probabilities calculated from the 256 x 256 matrix 11 of the original chain 
Xfc, which we denote by , are represented in blue. At time 0, the network 



.(0) 



.(0) 



1. The 



was put in state 1 (state 00000000) which is in Bl, so that 
state probabilities i-*^^ of the reduced chain Y^, calculated from the 4x4 matrix 
n, are shown in green. It is seen in Fig.[T4k that for p ~ 0.1 approximation Y^ 
is not good, in particular for basin B4. Results for p = 0.01 are presented in 
Fig. [T4b where it can be seen that the blue and green stairstep plots are almost 
identical, i.e. for p = 0.01, Y^ is a faithful coarse-grained representation of the 
NBN. 

Also compared the stationary probabilities calculated from X/j and Y/c (see 
1.4). For Xfe, we found (basins 1, 2, 3 and 4) 0.2856, 0.4577, 0.1259 and 0.1308 
when p = 0.1;0.2952, 0.4455, 0.1122 and 0.1472 when p = 0.01. For Y^, the 
stationary probabilites are independent of p (see below) and equal to 0.2963, 
0.4444 0.1111 and 0.1481. The maximum of the relative error is 13.26% when 
p = 0.1 and 0.94% when p = 0.01. Thus in the long run, the most populated 
basin is the one with the greatest size (120) and the smaller a probability (3/8). 

More generally, if chain Y^ is irreducible and aperiodic then its stationary 
state probability vector satisfies: 



zn 



1. 



(42) 



Rearranging the first equation in l(42|) we find: 



where 



Aai = 0, Zi + Z2 + ■ 



( -ai 

"12 
"13 



"21 "31 
-"2 "32 
"23 -"3 



\ "1_R "2fl "3i? 



"i?l \ 
"i?,2 
"i?.3 

-OLR j 



(43) 



From l(39|) . matrix is singular. The stationary state probability vector of Y^ 
is thus an eigenvector of A^, and the corresponding eigenvalue is 0. 

In addition, the stationary state probability vector can be obtained by cal- 
culating the limit 



^■'Notice that there is no a = basin (on > Vi). If on was null for some i, we would have 
everywhere in row i of 11 except at position i where we would have 1. Thus the reduced 
chain would be absorbing. 
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(a) 
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Time (fc) 




(b) 



50 100 150 200 



Time (fc) 

Figure 14: Reduction of a (8,2) noisy Boolean network. In blue: basin occupation 
probabilities z'^'^^ calculated from matrix 11 when initially the network is in state 
00000000 e Bl (state 1 of Xfe); in green: state probabilities ^ of Yfc calculated 
from matrix n (initial conditions: the chain is in state 1 of Y^). Transition probabil- 
ities TTij are estimated from (|40l). (a) p=0.1 (for the sake of clarity, the probabilities 
were interpolated linearly); (b) p=0.01 (stairstep plots and points omitted). 
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lim(/ + ^,)^ (44) 

A: — 'oo 

where each row of the Umit matrix is equal to the transposed stationary state 
probabiHty vector. 

5.1.2 Case p = 0.5 

When p — 0.5, the reduction is "exact". From (|36|) . the probabiHty to leave a 
basin of size i? is 1 — B /E. The probability to remain in such a basin is thus 
B/E. Since i/'s is geometric, one has: 

B- 

^iJ = (1 - if « ^ J' (45) 



E 



with 



T.Bi 



and BijE otherwise. Hence: 



B, 



l,2,...,i?. (46) 



The rows of 11 are thus equal. In fact, when p = 0.5, the stationary state 
is reached in one step whatever the initial conditions so that each row of 11 
gives the stationary state probabilities for the basins of the network. Note that 
since (5* — > as p ^ 0.5, these stationary state probabilities can be used to 
approximate the stationary state probabilities in a neighborhood of p = 0.5. 
From n, the mean sojourn time in basin i is expressed as: 

1 E 



E-B,, 



5.2 Continuous-time reduction 

The preceding sections deal with discrete-time homogeneous Markov chains. For 
such chains, the sojourn time spent in any state is geometric and thus memo- 
ryless. The continuous analog of the geometric distribution is the exponential 
distribution, which is also memoryless. Making the passage from geometric to 
exponential distribution leads to continuous-time homogeneous Markov chains. 
As we shall discuss, in the continuous-time representation of Markov processes, 
the state probabilities satisfy a system of linear ordinary differential equations. 

A continuous-time Markov chain {Yt, t > 0} is homogeneous if the transition 
probability from state i to state j in time interval {t, t + At) depends only on 
the length A t of the interval: nij{t,t + At) = Pr{Yt+At = j\Yt = i} =TTij{At). 
In this case (|Kleinrockl . [l975l) : 
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TTijlAt) = Uij At + o{At) if i j and 

l^nu{At) ^u,,At + o{At), (47) 
with Uij > the rate of transition from state i to state j ^ i and 



(48) 



The time spent in state i is exponentially distributed with parameter ua. 

Now if Yt is viewed as the continuous-time coarse-grained representation of 
an i?-basin NBN then the s tate p robabilities Zi{t) of Y* satisfy the following 
Master equation I Kleinrocl3 . 19751 ): 



■^Zi{t) = -UiiZ^{t) +'^UjiZj{t), z = 1,2, . . . (49) 

As stated in Proposition 3, to allow the passage from discrete to continuous 
representation, only one needs p to be sufficiently small. In other words, for 
sufficiently small p, one may use Yt instead of Y^ . A small p implies that the 
mean specific path c? = r is large compared to 1. Since ^J, > t, a small p also 
implies that ^ is large compared to 1. As an example, compare Figs. fT4k and 
fT4b . In the first case, d = 1.7558 and yU « 3 whatever the basin, while in the 
second case, d = 12.9441 and ^ is between 20 and 30 depending on the basin. 

Example 4 Let us illustrate the passage from discrete-time chain Xfe to continuous- 
time reduced chain Yt with the state diagram of Fig. [T] The reduced chain in 
this example has only two states so that from l|48p we get: un = ui2 = ui and 
U22 = W21 = U2. The equations for the reduced chain are thus: 

dzi 

—— = ~UiZi+U2Z2 

at 

= UlZl — U2Z2. (50) 

The expressions for the transition rates are found from l[40|l and l[47ll : taking At = 1, 
we get Ml — 7ri2/At = npai2 = np and U2 = 1121/ At = npa2i = np/2. 

Figs. [TSk and ITSb show Bl occupation probability versus time when p = 0.02 
and p — 0.002 respectively. The blue points represent the solutions of equation ((S} 
when initially the states of Bl are equiprobable and B2 is empty (for the sake of 
clarity, the points were interpolated linearly) while the green continuous curves are 
the solutions of system l[50|l when 21 (0) = 1 and 22(0) = 0. For p — 0.002 {fi = 
125.3756), the probability calculated from 11 decreases by small amounts and seems to 
vary continuously with time (see the enlarged portion in Fig. [TSbl so that continuous- 
time chain Yt may be used instead of Yt. Thus for sufficiently small p, system of 
differential equations l|50p may be used as a coarse-grained representation of the NBN. 

To end this example, let us estimate the relative error between the inverse of the 
mean sojourn time p,x = l//i and its approximation nap. The mean sojourn time 
for both basins and both p values are given in Table [21 For p = 0.02 one finds 
1/fii = 0.0749 and 1/^2 = 0.0369, to be compared to np = 0.08 and np/2 = 0.04, 
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which gives relative errors of 6.81 and 8.40%. For p = 0.002 one gets l//ii = 0.0079 
and 1/^2 = 0.0040, to be compared to np = 0.008 and np/2 = 0.004, which gives 
relative errors of 0.66 and 0%. 

6 Statistical fluctuations in NBNs 

Consider N replicas of a given BN, i.e. N cells expressing the same biochemical 
network, and suppose that each node of each replica may be perturbed with 
probability < p < 1 independently of time, of the other nodes and of the other 
replica0. 

The Markov chain model is a probabilistic model. Knowing the current 
state of a repHca, this model allows to compute the probability of finding the 
repHca in any given state of the network at any subsequent time. Therefore, 
even if the initial state of the replica is known with certainty, its trajectory in 
the state space of the network cannot be predicted with certainty. The same 
applies to the prediction of the number of replicas in each state of the network 
at any time. 

In order to illustrate the random behaviour of an ensemble of repHcas, ran- 
dom trajectories were simulated in the state space of the network of Fig. [T] by 
the Monte Carlo method. At time 0, N replicas were put in state 0010 G A2 then 
the number of replicas in each basin of the network at discrete times k = 1,2,... 
computed. The relative number of repHcas in B2 versus time is shown in Fig. [161 
for p = 0.02 and two values of TV. Graph (a) corresponds to = 10^ while 
graph (b) to iV = lO''. Each blue stairstep plot results from Monte Carlo 
simulations (one simulation is one trajectory of one replica), while each red one 
represents the mean solution calculated from ([8]). As can be seen from the two 
graphs, the uncertainty on the long-term behaviour of the ensemble is quite low 
in both cases {coefficient of variation: « 2% when A^ = lO'^ and « 0.7% when 
N = lO''). 

Remark 8 It is assumed that the total number of cells is conserved (cells do not 
proliferate and they are not lost): A^}*"' + Aj*"' = N \/k — 0,1, . . . 

Also calculated was the probability distribution of the time spent in basin 
82. The relative frequencies obtained from the Monte Carlo method are shown 
in blue in Fig. [iTlfor A^ = 10^^ and A^ = lO'' cases. The red points represent the 
exact frequencies calculated from matrix Q as explained in section 2 (Markov 
method) . For the sake of clarity, frequencies were interpolated linearly. Suppose 
that each Monte Carlo distribution in Fig. [T7| results from the measurement of 
A^ individual sojourn times. What would be the uncertainty on the mean time 
spent in B2 for each ensemble of cells ? The 95% confidence interval would 
be 27.2760 ± 5.8877% with the A^ = 10^ case and 26.9015 ± 1.9153% with the 
A^ = 10"* one. The exact mean sojourn time is 26.8067 (Markov method, see 
Table E]). 

Statistical independence between the replicas means that whether at time (fc + l) a replica 
has been perturbed or not does not depend on whether between and k other replicas have 
been perturbed or not. 
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Figure 15: Passage from discrete-time chain Xfc to continuous-time reduced chain 
Yt. Illustration with the state diagram of Fig. [TJ Ordinate: Bl occupation prob- 
ability; abscissa: time. Blue points: solution of ([8| when initially the states of Bl 
are equiprobable (the states of B2 are empty). For the sake of clarity, the relative 
frequencies were interpolated linearly; solid green line: solution of equations ([50]) 
with Ml = npai2 = np, U2 = npa2i = np/2, zi(0) = 1 and .22(0) = 0. (a): 
p=0.02; (b): p=0.002. 
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Figure 16: Statistical fluctuations in noisy Boolean networks. One considers N 
replicas of the Boolean network shown in Fig. [T](iV cells expressing the same bio- 
chemical network). Initially, the TV replicas are in state 0010 e 82. The graphs 
show the relative number of replicas versus time when p = 0.02 and (a) N = 10^ 
or (b) N = 10**. Blue stairstep plot: Monte Carlo method. Red stairstep plot: 
solution of matrix equation ([8|. 
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Figure 17: Statistical fluctuations in noisy Boolean networks. Probability distri- 
bution of the time spent in basin B2 of Fig. [T] when p = 0.02 and initially the 
N replicas are in state 3 (state 0010). Blue points: Monte Carlo method with 
(a) = 10'^ or (b) N = 10^. Red points: Markov method, exact frequencies 
(/z — 26.81, i5* = 0.11%). For the sake of clarity, the relative frequencies were 
interpolated linearly. 
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When N is fixed and p increases, the ampHtude of the fiuctuations in tps 
decreases as more and more rephcas have a short sojourn time. By comparison 
of the graph of Fig. [T6k and the two graphs of Fig. [181 the amplitude of the fluc- 
tuations in the relative number of replicas is not very sensitive to There is, 
however, a striking difference between Fig. [TSk and Fig. [TSb . When p is small, 
the relative number of replicas in B2 remains above or below its stationary value 
for long periods, i.e. the width of the fluctuations increases when p decreases. 
This means that in the low p regime, the long-term behaviour of the network is 
characterized by slow transitions between two states: one that corresponds to 
an "overpopulated" attractor and the other to an "underpopulated" one. The 
probability of crossing the stationary value during one time step was found to 
be 0.3111 when p = 0.1 and 0.0494 when p = 0.002 (0.1613 when p = 0.02). 
The maximum number of time steps the attractor remains overpopulated was 
16 when p = 0.1 and 318 when p = 0.002 (54 when p = 0.02). Similar values 
were found for the underpopulated case. 

Remark 9 (1) Since the number of replicas is conserved, when an attractor is over- 
populated, the other is underpopulated and vice versa. (2) These slow transitions oc- 
curing in the low p regime cannot be deduced from the Markov chain model. 

7 Reduction of a NBN to a two-state Markov 
chain 

We addressed the problem of reducing a chain to a two-state homogeneous 
chain by aggregating basins of attraction. 

Only (8, 2) networks with i? > 4 were studied. The aggregation process con- 
sisted in the following. When R was pair, R/2 basins were picked at random and 
aggregated, while when R was odd (i? — l)/2 basins were aggregated randomly. 
In both cases the remaining basins were aggregated, constituting the second 
state of the two-state chain. For each aggregated state then, we calculated the 
mean sojourn time /i with both types of initial conditions (the uniform and the 
random type) as well as the maximum deviation 5*. Results indicate that as 
p — > 0, /i does not tend to be independent of initial conditions neither does 6* 
tend to 0. Notice, however, that the reduction of Xfc worked well in some cases 
(some networks with some basin aggregations). 

8 Conclusion 

The reduction method for NBNs presented in this paper raises the impor- 
tant question whether biochemical networks can be reduced to (approximating) 
coarse-grained networks functionally equivalent to the original ones. Reduc- 
ing the complexity of biochemical networks could help in the analysis of cell 

^^Neither is the long-term B2 occupation probability: 0.6672, 0.6709 and 0.6711 when 
p = 0.002, 0.02 and 0.1 respectively. 
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Figure 18: Statistical fluctuations in noisy Boolean networks. Idem Fig.fTBb except 
that (a) p = 0.1 and (b) p = 0.002. In the low p regime, cells are trapped by 
attractors for a substantial time. Each attractor is alternatively overpopulated and 
underpopulated with regard to the stationary mean cell number. 
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responses to inputs (or cell fates) by neglecting molecular interactions while 
focusing on the higher-level processes that emerge from those interactions. A 
formally equivalent and very useful reduction theorem exists in electrical circuit 
theory which is Thevenin's theorem. 

9 Note 

This work is part of a manuscript entitled "Mathematical modeling of cellu- 
lar processes: from molecular networks to epithelial structures" written by F. 
Fourre. The complete manuscript contains five chapters. The first chapter is 
devoted to NBNs. The aim of the project is to propose a physical framework 
for describing cellular processes. Since 1st December 2008, F. Fourre has been 
working on a PhD thesis that is funded by the University of Luxembourg and 
supervised by Prof. Thomas Sauter. The thesis deals with qualitative modehng 
of signaHng networks. 

D. Baurain is a Postdoctoral Researcher of the FNRS. 
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